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I. INTRODUCTION 

The statistical mechanics of systems out of equilibrium is a rapidly evolving subject, due to the intensive research 
in the slow relaxation phenomena arising in several different contexts, such as coarsening systems, glassy and granular 
materials, colloidal systems etc. Understanding the basic mechanism underlying the slow relaxation is an issue of 
major importance. In particular, a key question is whether the large time scales are due to cooperative effects on 
large length scales. For non disordered coarsening systems this is certainly the case, since the observed power law 
relaxation can be directly related to the growth of the time dependent correlation length, or domain size [ij. In the 
case of disordered or glassy systems the establishment of such a connection is much more problematic, due to the 
difficulty of pinpointing the observables fit to the task. 

The use of the nonlinear susceptibilities has been advocated 0, |3, 0] as experimental or numerical probe apt to 
detect the heterogeneous character of the glassy relaxation and, possibly, to uncover the existence of the growing length 
scale responsible of the slowing down of the relaxation. However, this requires to establish clearly the relationship 
between non linear susceptibilities and correlation functions, in order to make sure what actually do the nonlinear 
susceptibilities probe. In other words, the problem of the derivation of the nonlinear fluctuation dissipation relations 
(FDR) in the out of equilibrium regime needs to be addressed. As a matter of fact, this has been one of the most fruitful 
lines of investigation in the field of slow relaxation, although mostly limited to the domain of linear response @, Q. 

In this paper we approach the problem on fairly general grounds. Within the framework of Markovian stochastic 
evolution, we bring to the fore the structural elements which are common to discrete and continous spins. We develop 
the formal apparatus necessary for a unified derivation of the FDR in the two cases and to arbitrary order. We also 
show, expanding on the work of Semerjian et al. Q, how the nonlinear FDR of arbitrary order can be derived from a 
fluctuation principle Q also in the off equilibrium regime. This allows to regard the FDR as a manifestation of the 
constraint imposed on the dynamics by the requirement of microscopic reversibility. 

The immediate application of the FDR is in the development of algorithms for the computation of the response 
functions without the imposistion of an external perturbation, the so called zero field algorithms. The numerical 
advantages of a zero field algorithm are remarkable. These have been illustrated and discussed in detail, in the linear 
case, in a recent paper |9i]. In the present paper we apply the zero field algorithm to the computation of the nonlinear 
response functions. We consider two cases, where knowledge of a nonlinear response function is required. The first 
one arises in the search for a growing length scale in the context of glassy systems, as mentioned above. The presence 
of quenched or self-induced disorder makes the linear response function short ranged, compelling to resort to the 
nonlinear ones. The second one deals with the extension of the effective temperature concept 10] to nonlinear order. 
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This is an important and difficult problem. Here, we consider it in the context of non disordered coarsening systems, 
showing the consistency of the effective temperatures derived from the linear and the second order FDR. 

The paper is organised as follows: the formal developments concerning the time evolution, the response functions, 
the FDR and the fluctuation principle are presented in sections 2, 3, 4 and 5, respectively. In section 6 the problem 
of the detection of a cooperative length trough a nonlinear susceptibility is addressed, while the effective temperature 
to nonlinear order is treated in section 7. Concluding remarks are presented in section 8. 

II. FORMALISM AND TIME EVOLUTION 

Let us consider a system in contact with a thermal reservoir, whose microscopic states are the configurations 
a — [ai] of the N degrees of freedom, discrete or continous, placed on the discrete set of sites « = 1, ..,N. Assuming 
a Markovian stochastic dynamics, the time evolution of the system is fully specified once the probability distribution 
P((T, to) at some initial time is given, together with the transition probability P{a,t\a' ,t') for any pair of times 
to ^ t' < t. Observables, or random variables, are functions A{a) defined over the phase space of the microscopic 
states, whose expectations are given by 

{Ait))=J2Aia)Pia,t) (1) 

a 

where P(a, t) = Y.^, P{a, to)P{(T' , to). 

In order to keep the notation compact, it is convenient to switch to the operator formalism [llj, whereby the 
microscopic states introduced above are the set of labels of the basis vectors 

k)=(8)l^»)' i = ^.-N (2) 

of a vector space. The single site states obey the orthonormality relation {(Ti\aj) = 5cri,ij'.5i,j- With this notation, the 
probability distributions P{a,t) become the time dependent vectors \P{t)) such that 

Pia,t) = {a\P{t)) (3) 

and the transition probability is associated to the propagator of the process P{t\t'), which is the operator whose 
matrix elements are given by 

P{<j,t\a',t') = {a\Pm\a'). (4) 

For stochastically continous processes, through the first order expansion 

P{t + At\t) = I + W{t)At + 0{At^) (5) 

where / is the identity operator, there remains defined the generator of the process W{t). Then, the pair \P{to)), W{t) 
contains all the information on the process, since the propagator can be written as 

P{t\t') =t( e/t' ''''^^A (6) 



where T is the time ordering operator. In differential form this is equivalent to the equations 

^^pm = wit)pm (7) 
-^pm = ^Pmwit') (8) 

from the first one of which follows the equation of motion of the state vector 

(9) 

Notice that the normalization of probabilities imposes on the generator the condition 

{-\W{t)^0 (10) 
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where (— | = X]<t('''I called the flat vector. We shall assume that at each instant of time the generator satisfies the 
detailed balance condition 



e' 



mt)^,(^t)^-mt) ^ w\t) (11) 



where 'H{t) denotes the Hamiltonian of the system, with a possible time dependence. This implies that the istantaneous 
Gibbs state 

\P('))0 = ^f~''^^'^\-) (12) 

is an invariant state in the sense that 

W{t)\P{t))f,^0. (13) 

The time dependent partition fimction is given by Z{t) — (— je"^^*-*-' |— ), as it follows from the normalization of the 
state {—\P{t))/3 — 1. From now on we shall adopt the notation |-)^ for the Gibbs states. 

Notice that Eq. ([T]) requires that observables correspond to diagonal operators. Each random function A{(t) is 
mapped into the operator A = A{a), where the operators ai are defined by 

a,\a)=a,\<7) (14) 

in such a way that the expectation ([T]) can be written as 

{Ait)) = {-\A\P{t)). (15) 

If W and \P)p are time independent, one can show that the multi-time expectations of observables in the stationary 
state obey the Onsager relation [l^]- Namely, if Ai, A2, An is a set of observables and i„ > in-i--- > ti an ordered 
sequence of instants of time, then 

(i„(i„)i„_l(t„_l)...i2(t2)il(il))/3 = 

(il(t„)i2(i« - it2 - h))...A„^i{ti + (t„ - <„_l))l„(tl))/3 (16) 

where 

(i„(t„)i„_l(t„_l)...i2(i2)A(il))/3 = 

{-\Ar,Pitn\tn-l)An-l-A2Pit2\tl)Ai\P)p. (17) 

Finally, using Eqs. ([7]) and ([5]), it is straightforward to show that the time derivative of a multi-time expectation is 
given by 

— {A4tn)..Ak{tk)..Ai{ti)) = 
Otk 

{An{tn).\Au,W{tk)]{tk).Mti)) (18) 

namely, the time derivative in front of the expectation amounts to the insertion, at the time tk, of the commutator 
[Afc, W{tk)] inside the expectation. This applies for all fc = 1, In particular, in the case of a single observable 

|(i(t)) = {-\AW{t)\P{t)) ^ {~\[A,wmP{t)) (19) 

where the second equality is a consequence of Eq. ([TU)) . It should be noticed that, in general, [A, W^(i)] is not an 
observable. 



III. RESPONSE FUNCTIONS 

Let us assume that the generator W{t) depends on time through an external field hi{t). Then, upon varying hi{t), 
there remains defined a family of stochastic processes and one may ask whether these stochastic processes are related 
one to the other. In particular, one would like to know whether the process with the generic hi{t) can be reconstructed 
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from the properties of the unperturbed process, the one with the particular choice hi{t) = 0. In order to answer the 
question, let us start from the statement that all the information in the process, with a given hi{t), is contained in 
the full hierarchy of the time dependent moments 



M,,,...,,„(ti,...,t„,[/i,(t')]) = {a,,{ti)...a^Jtn)) 
each of which is a functional of hi(t). Assuming analiticity, one can write the formal expansion 

M,,„,^{ti, ..,tn,[h,{t')]) = M,,„,^{ti, ..,tn) + 



(20) 
(21) 



oo 

V — V 

^ ml ^ 

m=l jl--jm 



dt[. 



dt' R 



rn ^-H.-inyi-.J, 



where il/i^. ^^ (ti, ..,i„) is the unperturbed moment, [t^^t) is the time interval over which the action of the external 
field is considered, and 



R 



{tl, .., tn', ti, .., t,„) 



<5'"M,,..,„(ti,..,t„,[/i,(0]) 



(22) 



h=0 



is the TO-th order response function of the n-th moment. Therefore, the question asked above can be positively 
answered if the response functions can be expressed in terms of quantities computable in the unperturbed process, 
namely if the FDR of arbitrary order can be obtained. 

Without loss of generality, we limit ourselves to work out the response functions for the first moment. For the 
responses of the higher moments there is nothing conceptually different, just the formalism gets more involved. As a 
matter of fact, in section IVll we shall deal with the second order response of the second moment. Coming back to the 
first moment, from 



M,(<, [h,{t')]) = {-\a,Phit\t^)\P{t^)) 
where Ph{t\t.ui) is the propagator in the presence of the field, follows 

,.. . .... 5"^pH{tK) 



\P{tu,)) 



(23) 



(24) 



h=0 



which obviously vanishes if tj ^ (i^,, t) for any j — 1, .., to. Let us write explicitcly the first two derivatives 



5hj^{ti) 



Jl dsw(s) dW{ti) 



and 



Shj,{ti)Shj^{t2) 



= T- 



f dsW(s) 

3 J till 



dW{ti) dW{t2) d'^Wjti) 
dh.Ati) dh,At2) ^ dhl{ti) 



(25) 



5(12) 



- P(f\f \^:^^p(f \f 1 ^^(^"0 p (f \f 



dhj„^itm) 



(26) 



where tM = max(tj), tm — niin(tj), jM,jm are the sites where the field acts at the times Im or tm, respectively, and 
5{np) — Sjn.jpS{tn ~ tp). The third order derivative is written down in Appendix I. 

In order to go further, we must specify how W{t) depends on the field hi{t) and that is where the distinction 
between discrete spins and continous spins comes in. In the following, in order to keep the derivation as simple as 
possible we shall specialize to the case of Ising spins with single spin- flip dynamics. The generalization to g-state spins 
(e.g. clock models) and to dinamical rules with conservation laws, such as Kawasaki spin-exchange, turns out to be 
straightforward. 
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A. Ising spins 

For Ising spins the state vectors jci — ±1) are represented by the column vectors 

k. = l)= (J), h: = -l)= (^J 
A generator of single spin flip dynamics is of the type 

I ^ 

W=-Y^W. (27) 

i=l 

where Wi has non vanishing matrix elements only between states which differ for the value of ai. The form of Wi is 
obtained imposing the detailed balance condition (jlip . Singling out the site i, the Hamiltonian can be written in the 
form 

U{t) = n{W\, t) + {h'^iia'W - h,it)}a^, (28) 

where af is the z Pauli matrix 

■ 1 

-1 

[a^ji stands for the set of all spins except for the one on the i-th site and hf^ = h^{[(j^]i) is the Weiss field on the 
site i. Inserting into the detailed balance condition (fTTj) . one finds the generator of the Glauber type 

W,{t) = (af - /)e'5[''^-'''(*)l*? (29) 



where af is the x Pauli matrix 



Hence, the derivatives are given by 



with 



1 

1 



a"VK(t 



Y = (-/3)'W,,(ti)(a|J" (30) 



(<j|J" = (^'/°^ - (31) 
For later reference, let us write here the identity 

obtained by adding the commutator and the anticommutator. It is straightforward to verify that the anticommutator 
is a diagonal operator and, therefore, is an observable which we will denote by 

Mt) = {i^-,W{t)}. (33) 

Due to Eq. (fTU]) , the average of the left hand side of Eq. ([5^ vanishes and Eq. for af can be rewritten as 

|(af(t)) = (B,(t)). (34) 



6 



B. Continuous spins 

In order to avoid confusion, here we shall denote by = [tpi] the set of the N continuous variables. Let us assume 
an equation of motion of the Langevin type 

^^^,{t)^BS)+V^{t) (35) 
where the drift Bt is related to the Hamiltonian, or free energy functional 7i [</?], by 

= 36 

Olfi 

rji (t) is a gaussian white noise with expectations 

(77.(0) = 0, {m{t)vj{t')) = 2T6,^,5{t - t') (37) 

and T is the temperature of the thermal reservoir. The generator of the corresponding Markov process is the Fokker- 
Planck operator 

W^'' = J2wr (38) 

i 

where 

Wr = -{Tp^^ + iB,[^]p, + D,m, (39) 

Bi [0] is defined through Eq. ((36|) and 



DM] = - 



obey the commutation relation 



and satisfy the equalities 



(40) 



The conjugated operators (pi and pi, defined by 

{ip'\if,\if) ^ ip,{ip'\if) (41) 



{^'mv) = (42) 

[if,,Pj]=i5uj (43) 
= (44) 



{-\Wr^Q- (45) 
The external field enters the free energy functional linearly 

=HM -^/i,(i)V5. (46) 

i 

yielding 

BhAvA^ B,[f] + Ht). (47) 
while -Di[(/?] remains unaltered. Therefore, the generator is changed into 

W^ = Wr-^h^m (48) 
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and the derivatives are given by 

^"^hT / —ipi, for n = 1 



dhf{t) \ 0, for n > 1. 
Obviously, the identity ([32| holds also in this case 



(49) 



and, using the definiton p9p together with the commutation relation (j43p . it is not difRcul to rewrite it in the more 
convenient form 
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^P^ = -{[W^,^,]+BM}. (51) 



Again, since the average of the left hand side vanishes, we get the analogous of Eq. ([M)) 

l{^,{t)) = {Bmt)) (52) 

showing that the anticommutator (|33p for Ising spins and the drift (|36p for continous spins play the same role in the 
evolution, thus justifying the choice of the same notation. 

IV. FLUCTUATION DISSIPATION RELATIONS 

Let us now return to the general treatment, valid for discrete and continous spins alike, with the operator a standing 
either for or for ip, depending on the context. Which is the case, will be specified whenever necessary. 

Comparing the left hand sides of Eqs. ([5^ and (|5ip with the first derivatives ([50]) and (HS)) . we can write the basic 
equation in this paper 



[w^,a^)] + B^[a]) (53) 



dW _ -/? 

where, in the case of continous spins, the superscript FP on the Fokker-Planck operator has been dropped and, as 
specified above, Bi stands for 

B,^ I {iS-*, W^}, for Ising spins 
* [ Bi[(p], for continous spins. 

We are now in the position to derive the FDR, by going through the following steps: 

1. according to Eq. (j24p . the response function is obtained by appropriate insertions of derivatives of the generator 
in between propagators, as exemplified in Eqs. (p5)l and (|26|) 



2. according to Eq. ([53]) . each first derivative of W can be repalced by the sum of the commutar and the drift 
operator Bi 

3. according to Eq. (jlSp . the insertion of the commutator amounts to a time derivative acting in front of the 
average. 

The above steps exaust all there is to do in the continous spin case, since there are no derivatives of the generator 
of order higher than the first. Higher derivatives do, instead, appear in the Ising case producing singular terms. In 
order to see how the procedure works in practice, let us carry out the computation up to second order. The explicit 
result for the third order response function is presented in Appendix I. 
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A. Linear FDR 

Let us begin with the simplest case of the hnear response function. From Eqs. (|24)) and ([25)1 follows 

dW{h) 



\cT,P{t\h) 



P{tl\t^)\p{t^)) 



h=0 



and, using Eqs. ([55)1 and p^ . this becomes the linear FDR 

d 



(55) 



(56) 



due to the appearence of unperturbed correlation functions of observables in the right hand side. This result has been 
known for some time for continous spins, see e.g. Ref. @], while for Ising spins has been derived for the first time in 
Ref. 0, where it has been exploited to develop the zero field algorithm for the computation of R'^^^^\t, ti), mentioned 
in the Introduction. 

If the system is in equilibrium, the averages in the right hand side are equilibrium averages and, invoking the 
Onsager relation ^TE\i . one gets 



{aKt)BM)p = {~\[a,,,W]P{t\h)a.,\P)i3 



d_ 
"dfi 



(57) 



after using space and time translation invariance in the last equality, Eqs. ([5^ or ([5T|l to eliminate Bj^, together with 
the normalization conditions (jlOp or Inserting this into Eq. (|56p . the equilibriimr fluctuation dissipation theorem 
(FDT) is recovered 



where we have introduced the pair correlation function 

Cy, (i, h) = Af,,, {t, ii) - M.(t)M,, ih) 



(58) 



(59) 



since it is clear that in equilibrium the one time quantities are time independent and the time derivative of Cij-^ (t, ti) 
coincides with that of Mij-^{t,ti). 



B. Second order FDR 



From Eqs. (|24p and (|26p . the second order, or two kicks nonlinear response function, is given by 

i?g,(i,tl,t2) = 



\a,P{t\tM) 



dW{t 



M) 



PitMK) 



dWiUn) 



h=0 

p{tM\p{t^))S{u) 



Ph{t^K)\P{Uo)) 



h=0 



h=0 



Substituting, next, Eq. (|53p for the first derivatives and the identity 

m^la^B^ + ^anatM 

for the manipulation of the singular term, eventually one finds the second order FDR 

4S.(^'*i'*2) = (/3/2)'{^^M.,.,,„(i,tM,i™) 
d 



(60) 



(61) 
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d 

d -I 

+ jr^—Mij^j^St^'tM ,tm) \6j^,,j,J{tM " tm) 

(62) 

where, it is worth recalhng, the last singular contribution in the braces is present only for Ising spins. 

Again, at stationarity the Onsager relation can be used to eliminate the -B's entering with the shortest time in favor 
of time derivatives. In so doing, the second and the fourth contribution in the right hand side become identical to the 
first and the third one, respectively. Furthermore, in the stationary state the last contribution in the right hand side 
vanishes, eventually yielding what we may call the second order FDT 

^Sl(*'^i'*2) = (/3V2)|^^C,,M,„(i,<M,i™) 

- ^{^:{t)BnAtM)<JjAtm))p'^ (63) 

where we have substituted the time derivative of the moment with that of the correlation function. Notice that there 
remains an ineliminable presence of B in the second term in the right hand side, which carries the information on the 
specific rule governing the time evolution. This is a distinctive feature of all FDR of order higher than linear, making 
them less universal than the linear one, which is B dependent only off equilibrium. 



V. FLUCTUATION PRINCIPLE 



In the following we show that the FDR obtained in the previous section arise as a consequence of a fluctuation 
principle Q. For convenience, the derivation is presented for the continous spin case, but it can be worked out 
along the same lines also for Ising spins. The idea of a derivation from the fluctuation principle was implemented by 
Semerjian et al. [3] in the stationary case. Here we extend the derivation to the more general off equilibrium case. 

Let us call experimental protocol the assigned time dependence of the external field h{t) in some time interval 
{tQ,tp). Then, from the detailed balance condition follows that the probability Pfs{[(p{t)]\(po, [h{t)]) of a path [<p(t)], 
taking place under the protocol [h{t)] and conditioned to the initial value ipito) = (pQ is related to the probability of 
the reverse path 

if{t)^ip{t), t^tp-t + to (64) 
under the reverse protocol [h{t)] and conditioned to ip{to) ~ ipp by 

P;3([(p(i)]|<^o,[Mi)])exp(-/3 r dth{t)^{t) 



= Ppmt)]WFAh{t)])e^^{mA^o)-nAvF)]} (65) 

where the subscript (3 is there to remind that the evolution takes place while the system is in contact with a single 
thermal reservoir at the inverse temperature /3. Multiplying both sides by ipi^FPi{ipo), where Pi{ipo) is an arbitrary 
probability distribution, and summing over the set C{to,tp) of all paths in the interval (to^tp) one finds 

[ d[p{t)]^,^pP0{Mt)]\(po,[hm(^^p{-l3 f dt h{t)^{t)\ Pi i^o) 

JC{to,tF) I Jto J 

= Zp.O f d[ip{t)]PI{cpo)e^'^°'^'^"^Pp{[m]\^FAHmV^,FP0,o{VF) 

(66) 
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where Pp^Q and Z^.o denote the equihbrium distribution and the corresponding partition function, in absence of the 
external field. Hence, the above result can be rewritten more compactly as 



(p,(tF)exp|-/3y dth{t)ip{t) 



l3,Q^I3,[h{t}] 



(67) 



where {•)j^0^[h{t)] stands for the average in the process starting with the initial condition P/, thereafter in contact with 
the thermal reservoir /3 and evolving with the protocol while {■) ^ {h(t)] ^^^^'^^ f'^'' ^he process in contact with 

the thermal resorvoir /3, starting with the unperturbed equilibrium distribution P^^o and evolving with the reverse 
protocol [h{t)]. 

The next step is to expand both sides in powers of [h(t)] about h{t) = and to compare terms of the same order. 
This is done in Appendix II, obtaining to zero order 



(68) 



where (•)/^^ q stands for the average in the off-equilibrium process starting with Pj and evolving in contact with the 
thermal reservoir in absence of the external field. At higher orders one gets 



S"{PiiipitF))e0^oi^('-)'>^,ito)) 



I3,0^f3,[hit)] 



Shj^{ti)...Shj^{tn) 



h=0 



P=0 •P(ii,..,i„„pb„_p+i,..j„) 



^^jn-p+i {tn-p+1 



h=0 



(69) 



where the second summation is over all the distinct permutations 

'Pi.hi --J jn-p\jn-p+i, ■■,jn) between the two sets of indeces (ji, j„_p) and {jn-p+i, ■■■,jn)- For instance, at the first 
two orders the above formula reads 



(5(P,(^(t;^))e^««(^(*-))^,(to)) 



/3,0-*/3,[?i(t)] 



h=0 



/3— (^.(t^)v.,,(ii))z^,,o+ jj^^ 



(70) 



/i=0 



and 



<52(p,(^(i^))e/5«o(^(*-))^,(to)) 



/3,0->/3,[h(t)] 



h=0 



d d 

p_d_ S{VtitF)(PjA'ti))i^fS,[h{t)] 
dti 5hj^{t2) 

S'^{ViitF))l~.f3,[h{t)] 



h=0 



^ d 6{ip^{tF)iP:j^{t2))l^l3.[h(t)\ 

dt2 Shj^{ti) 



h=0 



5hj,{ti )6hj^ (ta) 



(71) 



h=0 



From these two equations the off equilibrium FDR ((56)) and ((62)) can be recovered. Here we show this explicitely in 
the linear case, the extension to higher orders being straightforward. 

Recalling the definition ((^^ of the response functions and making the replacements tp ^ t and ~^ , Eq. ((70)) 
can be rewritten as 



d 



dh 



(72) 



h=0 
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According to Eq. ([49)) . the derivative with respect to hj^(ti) can be replaced by the insertion of —ipj^{ti) and, using 
Eq. ([5T|) . the second term in the right hand side can be rewritten as 



5hj,{ti) 



h=0 



P d 



^ ati \ / /3,o 



(73) 



Furthermore, since after setting to zero the external field the averages become equilibrium averages, using the Onsager 
relation we get 



f3,0^l3,lh(t)] 



h=0 



(3,0 



- Z^,o^(^.(i)i?,,(ti)Pz(^(i.))e''«°(^(*™»)^^^. 
The next step consists in the recognition that for an arbitrary observable A 



(74) 



(75) 



since the factor Z^,oP/(<<5(tu,))e''^°'^'^'-*™^' in the left hand side has the effect of undoing the equilibrium initial condition 
and replacing it with Pj. Hence, Eq. (|74p can be rewritten as 



^(3.0 

P_d_ 
2dh 



f3,0^f3,[h(t)] 



h=0 



(76) 



and inserting it into Eq. (|72p . eventually one finds 
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{ft{t)BjAtl))l^l3,Q 



(77) 



thus recovering Eq. (|56p . 



VI. NONLINEAR SUSCEPTIBILITY AND GROWING LENGTH SCALE 



In this section we consider the use of the response functions in the diagnostics of cooperative effects taking place 
over large length scales during the relaxation, referring to the case of Ising spin systems. A partial and preliminary 
account of the material in this section has been presented in Ref. [l3| • 

In non disordered coarsening systems, such as a ferromagnet quenched to or to below the critical point, the existence 
of a growing dynamical correlation length is well captured through the scaling properties of the equal time correlation 
function Cij{t) = Mij{t,t) [3], recalling that in this kind of processes Mi{t) = 0. In glassy systems, instead, quenched 
or self-induced disorder renders the two body correlation function short ranged, making it necessary to resort to higher 
order correlation functions. Attention has been particularly focussed on the four-point correlation function [isj 

^ii\'t^tw) — Miijj{t,t^,t,t^) — Mii{t,t^)Mjj{t,tj^) (78) 
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which describes the so called heterogeneities [l6j, namely the space fluctuations associated to the local time decor- 
relation. Although quite convenient in the numerical simulations, C^j\t,tu,) has the shortcoming of being hardly 
accessible in the experiments. Conversely, susceptibilities are more easily measurable and this has prompted the 
investigation of the nonlinear FDR. 

Here, we expand on the proposals 0, H, H, [13] to investigate dynamic scaling through higher order susceptibilities. 
Let us begin by considering the second order response of the second moment at equal times 



p(2.2) S^M,,{t,t, Mt']) 



Shj^{ti)Shj^ (^2) 



(79) 



Proceeding exactly as in the derivation of Eq. (I62p . the corresponding FDR is given by 

d d 

Rfimn^WiM) = {P/2f\ — — M,,,,,,^it,t,tM,tm) 

OlM Olm 

^ {at{t)aUt)al.^{tM)I3,M) 



dtM 

+ {(3^2) { (a,f (i)<T| {t)al {tM)B,^ (t^)) 
d \ 

(80) 

Looking at the above equation, it may seem farfetched to relate the properties of C'^^^ to those of R^^'^^ , since the fourth 
order moment appears explicitly only in the first term in the right hand side and in the singular term. Nonetheless, 
information on the existence of a growing correlation length can be obtained through a scaling argument, as it will 
be shown below. Let us consider the more manageable integrated response function 



where the Xi^'^'' the subtraction are the time integrals of the linear response function (j56p . i.e. 

X^''\t,t^)= I' dt,Rll^'\t,t,). (82) 

The reason for this subtraction and for the minus sign on the left hand side will be clear shortly. Assuming that 
eventually an equilibrium state is reached, from equilibrium statistical mechanics follows 

T\\imx^/\t,t^)=Cl^,g (83) 

where, for simplicity, we have taken (irf )eg = 0. The scaling relation C^.eg = ^■^^'^^''i^c.egd* — il/O; where ^(T) is 
the equilibrium correlation length, suggests a finite time scaling behavior of the form 



(84) 



where x = 4 — 2d — 2r] and z is the dynamical exponent. 

Another quantity, which has been recently considered in relation to C*-^-* , is the third order integrated response 
of the first moment 

[]'''\t,t^) = -- dhdt2dhR^^if]{t-tiMM) (85) 



XlJ V 7 - / 2 



where Rl-ijj is given by Eq. (|llip of Appendix I. Again, from the large time result 

,(1,3), 



T^^hm x,V'^^(i,t.) = C4,, (86) 
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we may infer the scahng behavior 

TH^'^^(M.) = rG(t^,^,^). (87) 

Notice that the prefactor in the definition (|85p has been arranged in such a way that the large time limits of T^x'^'^"* 
and T'^;^(i'3) are the same. 

These patterns of behavior have been checked numerically in the one-dimensional Ising model. The simulation 
has been carried out through standard Montecarlo techniques with Glauber transition rates, where Bi — af — 
tanh(/3 X]<j> ^v^j) and the sum runs over the nearest neighbours < j >i of i. Taking Jy = 1, we have prepared the 
system in a high temperature uncorrelated state and then quenched it to the final temperature T at time t — 0. In 
this case z = 2 [3, and C(T) = -1/ ln[tanh(l/r)] [l|. The integrated response functions xf/\t,0) and x!j'^^(i,0) 
have been computed using the FDR (|80lllip . following the zero field method of Ref. 0. It must be stressed that, 
due to the extremely noisy nature of the response functions, the numerical computation of these quantities through 
the FDR is by far more convenient than the computation based on the application of a small external field. In fact, 
in addition to the incomparably better signal to noise ratio, the — » limit is built in the FDR. 

In order to verify the scaling relations and (|87p . first of all we have taken ^ in order to reduce the variables 
from three to two in the right hand sides. Then, recalling that for the Id Ising model rj — I implies x ~ and varying 
the temperature and the distance in such a way to keep \i — j fixed, it is matter of showing that T^x|^'^^(t, 0) and 
r^Xij''^^(^i 0) are functions only of t^l^ j^. This is shown, with good accuracy, in Fig. ([T|) for the quenches to the three 
final temperatures Ti = 0.6, T2 = 0.7572 and = 1.0239 and with six different values of \i-3\li. Both T"^ ■,^) 
and r^x!j''^''(ii 0) grow from zero to the same asymptotic value C^^ — exp{— 2|i — on the same timescale, as 

expected. Note that the data for y;^^^'^^ are much more noisy than those for x*-^'^-*- Since the same numerical resources 
have been allocated in the computation of each of these two quantities, we conclude that investigations based on x^^'^'' 
are more efficient, at least numerically. 

We now discuss how x*^^'^'' can be effectively used for the measurement of a cooperative length in disordered 
systems. In order to improve the statistics further, it is convenient to consider the = component of the space 
Fourier transform 

x£o (*'^-) = (l/iV')ExS''^(^:i-) (88) 

which, using Eq. ([84]) . scales as 

In Ref. [ij] we have computed numerically this quantity, with = 0, in the Edwards- Anderson (EA) model with 
Hamiltonian 7i = ^j^- Jij&^aj and d — 1,2. The data have been analysed as follows: for large t the curves saturate 
to the equilibrium value Xg^Q^ ^ ^4-(i-2?7 Using the known values of 77, one can extract ^. In the off-equilibrium 

regime ^1/^ < C the growing correlation length L{t) ^ t^l^ is expected not to depend on ^. Enforcing this condition 
from Eq. ^ one must have J^(ii/V^,0) - (^/L(t))2''+''-4, which yields 

xg(t,0)^i(0^-'^-2^ (90) 

This allows to determine Lit). After doing this, we checked for the data collapse by plotting C '*^'^^^''Xfc=o^(^? 0) vs 
L{t)/^ for all the temperatures considered (see Figs. 1213]) . 

Let us consider first the d = 1 EA model with bimodal distribution of the coupling constants Jij = ±1. This system 
is considered in order to test the method, since it can be mapped onto the ferromagnetic system, just considered above. 
Moreover, in this simple case, in addition to X^^o (^'0)i '^'^'^ can check the scaling of the equal time structure factor 
Ck=o{t)i obtaining independent determinations of L{t) and ^ to compare with. This shows that the sets of data for L{t) 
and ^, obtained in both ways, are in agreement with each other and with the analytical behaviors up to the numerical 
uncertainty. The data collapse of Xfc=o^(^iO) and of Ck^oit) are shown in Fig. [21 Here, one clearly observes the off 
equilibrium regime, characterized by the power-law behavior of X^^o ^k^o (t) with exponents A — d — 2r] and 

2 — ?7, respectivly. In the large time regime equilibration takes place with the convergence of Xk^o ('^^ ^) Ck^Q{t) 

f2 2) 

to Xk=0,eq and to Ck=0,eq{t). 
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FIG. 1: (Color online). Upper panel: The quantity T^Xij'^' 0) is plotted against t^^^/£, [z = 2) for six fixed values of 
r = {r = 0,0.2855,0.5711,0.8566, 1.1421, 1.4277 from top to bottom) and three temperatures Ti = 0.6, r2 = 0.7572, Ts = 

1.0239 such that ^{Ti)/£_{T2) = 2 and ^{Ti)/^{T3) = 4. Lower panel: Same for r^Xi!'^^*, 0). 
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FIG. 2: (Color online). Data collapse of Xfe=o^ {Ck^o in the inset) for several temperatures in the d - 
The dashed lines are the expected power-laws in the non equilibrium regime. 



1 EA model (77 = l,z = 2). 



After this explicit verification, we have turned to the c? = 2 case, where the independent information on the structure 

(2 2) 

factor is not available. Both with bimodal and Gaussian distributions of J^, , the behavior of £ extracted from y- , 

k=0,eq 

using rj ^ 17 1, has been found consistent with previous results [17lll8l|. The non-equilibrium behavior is compatible 
with a power law L{t) ~ with a temperature dependent exponent in agreement with z{T) ~ 4/r, as reported 

in Ref. [l^. The data collapse of Xk=Q ^) shown in Fig. [3l Notice also the additional collapse of the curves with 
bimodal and Gaussian bond distribution, further suggesting that the two models may share the same universality 
class at finite temperatures [l7|. 



VII. EFFECTIVE TEMPERATURE 



One of the most interesting developments in the study of the linear FDR out of equilibrium has been the introduction 
of the concept of effective temperature T^f f (lo| . The idea is that the off equilibrium behavior observed during slow 
relaxation can be accounted for by the separation of the time scales for different subsets of degrees of freedom. Each 
of these is regarded as in equilibrium with a different virtual thermostat at some appropriate effective temperature, 
which depends on the time scale and is different from the physical temperature of the real thermostat driving the 
relaxation. The value of T^ff can be inferred by forcing the off equilibrium linear FDR in the form of the equilibrium 
FDT. Although appea ling , this idea has turned out not to be applicable tout court, since T^ff might turn out to be 
observable dependent [20|. Nonetheless, with the proper caveats, the concept remains quite useful and suggestive. In 
this section we make a preliminary exploration of another open end in the important question of how general T^f / can 
be, investigating whether it is possible to extend to the nonlinear FDR the effective temperature concept, consistently 
with what it is done in the linear case. 

Let us first recall how T^ff is defined from the linear FDR. For definiteness, the Ising spin case will be considered. 
Writing explicitely the time integral in Eg. ([5^ . one has 



,(14) 



(t, t^) — — 



dti 



(91) 



Assuming Mi{t) = throughout the dynamical evolution, and replacing Mii{t,ti) with the autocorrelation function 
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FIG. 3: (Color online). Data collapse of xi^o' 
(filled symbols) bond distribution, with z{T) — 
regime. 



for several T in the d — 2 EA model with bimodal (open symbols) or Gaussian 
4/r and rj = 0. The dashed line is the expected power-law in the non equilibrium 



C{t,ti), the quantity 



* d 



^^^\t,t^)^ I dti—Cit,h) = l-C{t,t^) (92) 



for fixed is a monotonously increasing function of time, which allows to reparametrize t in terms of i/'^^-' and to 
write Xi^'^^(i, ^m) in the form 

X^''\t,t^)^X^''\^^'\t^). (93) 

In equilibrium, where time translation invariance holds, the dependence on disappears and the parametric repre- 
sentation becomes linear 

xr^'^(V'(i))=/3V(i) (94) 

with the obvious consequence 

^-^V^w • ^^^^ 

Off equilibrium, the parametric representation won't be linear and an effective temperature can be defined by the 
generalization of the above relation 

with Peff = ^/Teff- 

In order to see how T^ff ^ T arises in a simple context, let us consider the relaxation to a low temperature phase 
characterized by ergodicity breaking and, therefore, by a non vanishing Edwards- Anderson order parameter qea- 
In particular, let us think of the already mentioned coarsening process, like in a ferromagnet quenched to below the 
critical point and relaxing via domain growth. In that case qea coincides with the spontaneous magnetization squared 
M2 . 



17 



As oo, the separation of time scales takes place. The short, or quasiequilibrium, time regime holds for 

C > A/g^, that is ■ip^^'> < 1 - M^^, while the large time scale sets in when C < M^^, or ■ip^^'> > 1 — M^^. The existence 
of this latter regime makes it evident the failure of equilibration, even in the t^^ ^ oo limit, since the auto corelat ion 
function falls below the Edwards- Anderson plateau. The behavior of ip^^\ obtained from numerical simulations of 
the Ising model in d = 2, is shown in the inset of Fig. [4l Starting from zero, there is a fast growth in the short time 
regime, followed by a plateau, more evident for large tw, and finally there is convergence, with the power law behavior 
1 — ijj^^\t,tw) ~ t~^/^, toward the asymptotic value V''-^'' = 1- Notice that A is the Fisher-Huse exponent [2l| and 
that the plateau flattens over the asymptotic value 1 — qea as t^t, — > oo. 

Correspondingly, the integrated response function can be written as the sum of two pieces (22| 

xf^'^ i^^'^ ,t.)- xi:-'^ {^^'^ ) + xiY^ {^^^ , t.) (97) 

where xit'^^ is the stationary contribution arising from the equilibrated bulk of domains, while Xag^"* is the aging 
contribution due to the off equilibrium domain walls. The stationary contribution obeys Eq. (|94p in the short time 
regime, saturates to its equilibrium value and remains constant in the large time regime, while the aging contribution 
vanishes as tm — > oo according to 

X^^g'\^^'\t^)=t--F{rj:(^^) (98) 
where a > [13] ■ Hence, the full response function obeys the asymptotic form 

lim v(^^^)(^(i) t)-{ ^ '^^'^ ^ ^ - . (99) 



which, on account of Eq. ([M]) . leads to 



V R nil) , \ r /3, for < 1/^(1) < 1 - Af2 

Namely, the effective temperature coincides with the physical temperature in the short time regime, where the sytem 
appears equilibrated, while it is drastically different from it in the off equilibrium large time regime. 

Let us now carry out the parallel analysis on the fc = component of the second order integrated response function 

43 (*'^-) = t^e/ ^^i/ dt,i^^j(t,t-t,,t,). (101) 



Notice that, although we use the same notation, this quantity differs from the one in Eq. ([88)1 . because of the overall 
sign and of the absence of the subtraction. Let us introduce the quantity 

+ (af(t)(7|(i)^,(iM)4(i™))] (102) 

with properties similar to those of ^^'^''{tjtw), as shown in Fig. |31 The main feature is the monotonous increase 
from zero to an asymptotic value well above the limiting value (1 — Qea)'^ that one would get from the equilibrium 
calculation. This is the value at which a plateau develops as gets large, signaling the separation of time scales. 
Using '0*^2^ to reparametrize the time t, from Eq. (|63p follows that at equilibrium the FDR becomes linear 

4=0 (103) 

Hence, by following the same reasoning as in the linear case, in the off equilibrium regime an effective temperature 
can be introduced by the analogue of Eq. ((96l) 



The question, now, is whether the two /3eff defined by Eqs. ([96]) and (|104p are consistent or not, that is whether the 
equality 

lim p,ff{^p^^\t^)^ lim /3e//(V'''\t») (105) 
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FIG. 4: (Color online). The quantity V'^* is plotted against t - i„ for t„ = 100, 500, 1000, 2000, 3000, 4000, 5000 from top to 
bottom. Data refer to a 2d Ising system of size 1800^ quenched to T = 2 (T^ ~ 2.269). The dashed horizontal line is the 
large-tu, height of the plateau {y = {1 — qsA)^)- In the inset the quantity tp^^^ is shown (same t^j of the main figure). 



holds or not. This is a difficult question to answer in general, we shall limit to the consideration of the particular 
coarsening process analysed above in the linear case. 

The short and the large time scales, in terms of ip'"'^\ correspond to tp^^'' smaller or larger than (1 — Af^^)^, 

respectively. Writing X^^^q eis the sum of two pieces, as in Eq. ((97)) . 

' = ) + ^S''^ i^^'^ ' (106) 

the same considerations made on xit'^'' apply exactly to x^^t''^\ since this is an equilibrium contribution. Namely, 
after obeying Eq. (|103p in the short time regime, saturates to the equilibrium value (1 — M^^)"^ and then remains 

constant in the large time regime. For x'af \^'^^\t^) there are no previous results to rely on. We have, then, 
measured xig'^"* {ip^'^\tw) numerically in the quench of a two dimensional Ising model below Tc evolving with Glauber 

(2 2).. 

dynamics. The agmg contribution Xag has been isolated using the method based on the no-bulk- flip dynamics 
discussed in [H [H [13 . 

In order to check if a scaling form of the type ([M)) is obeyed 

x^^fH^^'\t^)^t-^'m/t.) (107) 

we have plotted Xag {'^^'^w) foi' ^ fixed value of t/t.yj against tyj^ as shown in the inset of the upper panel of Fig. [S] 
From the observed power law behavior we have extracted the exponent a2, finding values in the range [0.59 — 0.62]. 
We have, then, carried out the data collapse by plotting f^Xag {t,tia) against t/tw For large the collapse (Fig. [5]) 
is quite good, confirming that the scaling form (|107p is obeyed with an exponent 02 — 0.61 — 0.62. For small values of 
tu) and t/tw deviations are observed due to preasymptotic effects, similarly to what was already observed in the linear 
case jl^ll^l- Notice, also, that for large t/t^ one has a power law decay of the scaling function F2{t/tyj) ^ {t/tw)~°'^ 
with the same exponent a2 entering Eq. p07p . exactly as it was observed in the linear case [27j . 






0,1 
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FIG. 5: (Color online). tl^Xag ii^,^^) is plotted against t/t^ for t„ = 100,200,500,1000,2000,3000,4000,5000 (from bottom 
to top) and three temperatures T = 1, 1.5, 2 (upper, central and lower panel). The exponent a2 is extracted as discussed in the 
text, finding 02 = 0.61 for T — 1, and 02 = 0.62 for T = 1.5 and T = 2. The dashed lines represent the asymptotic behavior 
Xag {t,tiu) ~ {t/tiu)~"'^ for large t/tnj- In the inset the quantity Xag {t,tw) is plotted against with t/t^ fixed. 
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FIG. 6: (Color online). The parametric plot of xi=o'' ('-i ''™) against /S^?/)^^' {Xi^'^\'tj'tw) against /St/;'^' in the inset) is shown 
for r = 2 and t„ = 100, 200, 500, 1000, 2000, 3000, 4000, 5000, from top to bottom. The dashed line is the expected asymptotic 
behavior. 



In conclusion, like in the linear case, the existence of the scaling behavior (|107p with a2 > implies that the aging 
contribution vanishes asymptotically, yielding the analogue of Eq. (|99p 

(2.2) , . r for ^(^) < (1 - M^^r 

,i™o.^£=o - [ f3^1 - M2^)2, for (1 - M^^r < i^^^y ^ ^ 

The approach to this asymptotic behavior is shown in Fig. [S] Hence, using the definition p04p . we find 

The comparison with Eq. (|100p suggests that the consistency condition (jlOSp is satisfied. 



VIII. CONCLUSIONS 



In this paper we have derived the off equilibrium FDR of arbitrary order for systems evolving with Markovian 
stochastic dynamics. The main effort has been to put the FDR in the same form for both continous and discrete 
spins. In order to stress the generality of the result, we have also shown how the whole hierarchy of FDR can be 
made to descend from the fluctuation principle. Once the FDR are available, response functions of arbitrary order 
are expressed in terms of unperturbed correlation functions of observables. The payoff is in the development of simple 
and efficient zero field algorithms for the numerical simulations. 

As an application, we have considered the problem of detecting the existence of a growing length in those cases, 
like in glassy systems, where standard methods based on two-point correlation functions and the corresponding linear 
response functions are of no use. In these cases the simplest object carrying useful information, in principle, would 
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be a four-point correlation function which, however, is not directly accessible to experiment. Instead, experimentally 
accessible are the nonlinear response functions involving the four-point correlation function through the nonlinear 
FDR. The choice of which response function and, therefore, of which FDR to use is not univocal, once the realm of 
the nonlinear response functions is entered. Then the choice is matter of convenience. We have made the proposal 
to use the second order response of a two-point correlation function, rather than the third order response of the 
magnetization, as advocated elsewhere in the literature. We have, then, demonstrated the numerical advantage of our 
choice through the implementation of the zero field algorithm. 

Finally, we have made a first step into the important but difhcult problem of the definition of the effective tem- 
perature through the nonlinear FDR. We have considered the domain coarsening process ensuing the quench of a 
ferromagnet below its critical point. Indeed, in that case we have found that it is possible to extract from the nonlin- 
ear FDR an effective temperature which is consistent with the effective temperature obtained from the much studied 
linear FDR. 



IX. APPENDIX I 



Proceeding like in the derivation of Eq. (j26|) , the third order derivative of the propagator is given by 

+ Ph (t I tM ) (^^^ i^I I ^™ ) Qf^. [trn \ tw ) d{tM ^ tj) 

+ Ph{t\tM) I Ph{tM\tl) ^Tff.^ Ph{tl\t^)5,,,jJ{ti - tra) 

+ Pkit\h)^^^PH{h\t^)6{12)5i23) (110) 

where tM — max(ij), t„i = min(tj), <tj <tM, jMijiijm are the sites where the field acts at the times i^/, tj or 
tm, respectively, and S{np) — 5j^,j 5{tn — tp). Inserting this into Eq. ([24]) and using Eq. ((53|) . we get the third order 
response of the first moment 



\ d 

-— — — (CT; {t)aj^, (tM) CTj, {ti)Bj^{tm)) ~ „, „, {cTt {t)(TjM (J^M ) Bj, {tj ) (Jj^^ (t„ ) ) 
OlMCltl aiMOtm 



+ -g^{^^it)BJ^AtM)^JJitI)BJ^itn^)) + {<^^{t)(T, M)B,At l)Bj^itm)) 

- {<t)BuAtM)Bj,{ti)B,^^{tm))] + ^[[-^{cr^{t)aJ,{tI)B,^{tI)a,^{U,)) 
d d ^ ^ . . 

-Mij„tj,j^{t,tM,tl.,tm) - {(7i{t)dj,{ti)Bjj{ti)Bj^{tm)) 



dt,n dti 

d " \ f ^ 

-Q^^{^i{t)^3M{tM)^ji{ti)Bj^^{Ua))j5j^aj5(tM - tl) + y-Q^^{odt)cFj^,{tM)cFj:{tl)Bj,{ti)) 

d d . - . - 

-MijMji],^ {t, tM,ti, tm) - {ai{t)Bj^j {iM)oj, (ti)Bjj (ti)) 



dtM dtr, 

■^{Mt)Bj,AtM)aj,iti)a,^ itm)))6j,j,J{ti ~ ^(Am,,, (t,ti) 

{a,it)B,,{t,)))6,,,Jj,j,diti - t2)5{t2 - t,). (Ill) 
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At stationarity this becomes 

/33 r 93 



(1 3) r r ^ 



^2 



dtMdt„i 

+ ^{^r it)Bj,r {tM)Bn {im))] 



Otra Oti 

d 



^ {t, tM,tl,tm)^6j^jjjS{tM — tl) 



(112) 



+ P''^^M,j,it,h)6j,jJj,jJih-t2)Sit2-h). 

It should be recalled that the singular terms in the last two equations are present only in the Ising spin case. 

X. APPENDIX II 

The expansion of the left hand side of Eq. ([S7)) can be done in two steps. Expanding first the exponential we get 

(p,(tF)exp(-/3 /" dt h{t)ip{t)\\ 

= X! i~ X! / d,ti...dtjn{fz{tF)>fjAtl)---'fjm{'tm))l^fJ,[h{t)] 

hj,{ti)...hj^{tm) (113) 

and expanding the individual averages 



^ — n i n, n -J to 



=0 ^ qi...qp 



h=0 



h,,{t[)...h,^{t'^) (114) 
all together these two contributions give 

^,{tF)expl-f3 f dth{t)<f{t) 



E E / dt,...dt^dt[...dt'^ 



'to J I I~,li\h{t)\ 



m=Op=0 ?l---9p 



h=0 



Ito " Sh,,{t[)...Sh,^{t'^) 

h,Ati)-h,^{Un)h,AA)-\AO- (115) 
Reorganizing the double sum 

EE rnlpl E E "EE (n-py.pl E E (11^) 
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the above result can be rewritten as 



ipiitp) exp < — /3 / dt h{t)ip{t) 



oo ^ n 



n. 

n—l p—0 



7.1 



E E 



SP{ip^{tF)iPjAh)---^j„^p{tn-p))l-.0,[h{t)] 



{ipi{tF)) 1^/3,0 

dti...dtn-pdt'i...dt'p 

hj,{tl)...hj,^_^{tn-p)hq,{t[)...hq^{t'p). 



h=0 



(117) 



Notice that the combinatorial factor in the square bracket gives the number of the distinct permutations among the 
two sets of indeces (ji, ■■■,jn~p) and (qi, q,p). 

Going over to the right hand side of Eq. (j67p and introducing the shorthand 



/3,0-/3J/i(t)] 



/3,0^/3,[/i(t)] 



(118) 



one gets 



hj,{ti)...hj,^{tn) 

and, since hj(t) — hj{t), this can be rewritten as 



«=l"-ji...j„-^*o 



Shj^{ti)...Shj^{t„) 



h=0 



(119) 



E^ E / " dt,...di„ 



/3,0^/3,[/i(t)] 



hj^{ti)...hj^(tn) 



5hj^{ti)...5hj^{tn) 



h=0 



(120) 



where (■)/3,o stands for the equilibrium average at the temperature (3 and without external field. Therefore, comparing 
with Eq. (fTT7)l . one arrives at Eqs. (UHl) and (ISg)) . 
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